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Abstract 

In  this  paper,  the  Keller-Box  method  has  been  applied  to  the  coupled  one-dimensional  electrochemical  transport  equations  in  order  to  simulate 
lead-acid  batteries  numerically.  The  advantages  and  disadvantages  of  this  method  have  been  discussed.  The  results  indicate  that  the  Keller-Box 
method  is  a  suitable  method  for  integration  of  electrochemical  transport  equations  both  in  integrated  and  multi-region  formulation.  The  boundary 
conditions  and  interface  conditions  (in  the  case  of  multi-region  approach)  can  be  implemented  easily  and  require  no  special  routine  (for  off-diagonal 
terms).  In  addition,  the  effect  of  acid  concentration  dependency  of  open  circuit  voltage  has  also  been  investigated. 
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1.  Introduction 

Computational  fluid  dynamics  (CFD)  can  be  used  to  solve 
the  governing  equations  of  battery  dynamics.  This  makes  sig¬ 
nificant  progress  toward  an  insightful  understanding  of  battery 
performance  and  detailed  characteristics  of  its  behavior.  To  have 
a  good  understanding  of  battery  performance,  different  design 
parameters  of  the  lead-acid  batteries  should  be  analyzed.  Tradi¬ 
tionally,  these  design  parameters  are  evaluated  experimentally 
which  is  time  consuming  and  costly.  Numerical  simulation  of 
the  battery  is  an  alternative  to  evaluate  the  battery  performance 
and  optimizing  the  design  parameters.  Another  advantage  of  the 
CFD  modeling  is  that  it  can  be  used  as  a  virtual  model  for  bat¬ 
tery  with  which  the  necessary  parameters  for  dynamic  modeling 
must  be  obtained.  In  this  method,  the  partial  differential  system 
of  equations  governing  the  behavior  of  the  battery  dynamics  is 
solved  numerically  using  advanced  CFD  techniques. 

The  governing  equations  of  the  battery  dynamics  [1]  have 
been  developed  in  different  forms.  Gu  et  al.  [6]  proposed  a 
multi-region  system  of  equations  to  simulate  the  battery  dynam¬ 
ics.  In  this  model,  each  region  was  studied  separately,  i.e.  for 
each  region  (Fig.  1)  a  system  of  transport  equations  was  given. 
To  relate  the  regions  at  their  common  boundaries,  a  set  of 
interface  conditions  was  also  proposed.  This  system  of  equa¬ 
tions  was  solved  numerically  with  finite-difference  method.  To 
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keep  second-order  accuracy  at  interfaces,  off-diagonal  elements 
appear  in  the  dicretized  system  of  equations  destroying  the  tridi¬ 
agonal  nature  of  the  resulting  matrix.  Fan  and  White  [2]  proposed 
an  auxiliary  routine  called  MB  AND  to  treat  the  off-diagonal  ele¬ 
ments  without  losing  the  second-order  accuracy. 

Later,  Gu  et  al.  [7]  introduced  an  integrated  formulation  for 
battery  dynamics.  In  this  approach,  the  whole  battery  was  con¬ 
sidered  as  a  model  volume  and  the  transport  equations  were 
derived  for  the  whole  cell  volume.  With  this  formulation,  the 
interface  condition  is  no  longer  necessary.  The  proposed  system 
of  equations  was  solved  numerically  by  means  of  finite- volume 
method. 

In  this  study,  the  Keller-Box  method  is  used  for  numer¬ 
ical  integration  of  transport  equations  of  battery  dynamics. 
This  method  can  be  implemented  on  both  approaches.  The 
Keller-Box  method  is  an  implicit  method  which  is  second- 
order  accurate  in  both  time  and  space.  Since  only  two  points  are 
involved  in  discretization,  a  nonuniform  grid  can  be  used  without 
any  difficulty.  In  addition,  the  Keller-Box  method  works  both 
with  the  unknown  functions  and  their  derivatives  at  each  grid 
point  simultaneously.  This  property  would  ease  numerical  for¬ 
mulation  as  well  as  implementation  of  boundary  conditions  and 
interface  conditions  (if  any).  The  disadvantage  of  the  method 
is  the  size  of  coefficient  matrix  increases  due  to  introducing  the 
derivatives  as  unknowns  which  increases  the  computational  cost. 
But,  once  the  solution  is  obtained,  the  derivatives  of  unknown 
functions  are  available  and  no  new  further  formulation  is  needed 
to  obtain  the  derivatives  from  the  unknown  functions  as  is  done 
in  the  other  numerical  methods.  The  results  are  compared  with 
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Nomenclature 

a  coefficient 

A  specific  electroactive  area  (cm2  cm-3) 

c  acid  concentration  (mol  cm-3) 

D  diffusion  coefficient  (cm2  s_1) 

F  Faraday  constant,  96,487  C  mol-1 

I  applied  current  density  (A  cm-2) 

j  transfer  current  density  (A  cm-2) 

k  conductivity  of  liquid  (S  cm- 1 ) 

m  molality  of  the  acid  (mol  kg- 1 ) 

t  time  (s) 

Greek  letters 
s  porosity 

a  conductivity  of  solid  matrix  (S  cm-1) 

0  electric  potential  (V) 

Subscripts  and  superscripts 
D  pertinent  to  diffusion 

eff  effective,  corrected  for  tortuosity 

ex  exponent  in  the  effective  property 

1  liquid  solution 

o  initial  value 

s  solid  phase 


the  results  of  previous  investigators.  Furthermore,  the  effect  of 
acid  concentration  dependency  of  open  circuit  potential  on  cell 
voltage  and  acid  consumption  is  studied. 

1.1.  Governing  equations 

A  typical  lead-acid  cell  is  shown  schematically  in  Fig.  1  and 
consists  of  the  following  regions:  a  lead-grid  collector  at  v  =  0 
which  is  at  the  center  of  the  positive  electrode;  a  positive  Pb02 
electrode;  electrolyte  reservoir;  a  porous  separator;  a  negative  Pb 
electrode;  finally,  a  lead-grid  collector  at  x-l  which  is  at  the  cen¬ 
ter  of  the  negative  electrode.  The  positive  and  negative  electrodes 


Region  1  Region  2  Region  3  Region  4 

Positive  Electrode  Reservoir  Separator  Negative  Electrode 


jc =0  x=I 


Table  1 

Governing  equations  of  the  battery  dynamics 

Conservation  of  charge  in  solid  V  •  ((jeffV0s)  —  Aj  =  0 

Conservation  of  charge  in  liquid  V  •  (&effV0i)  +  V  •  (&^ffV(ln  c))  +  Aj  =  0 

Species  conservation  =  V  •  (DeffVc)  +  ai 


consist  of  porous  solid  matrices  whose  pores  are  flooded  by  a 
binary  sulfuric  acid,  H2SO4.  The  model  is  assumed  to  be  one¬ 
dimensional  perpendicular  to  the  face  of  the  electrode.  During 
charge  and  discharge,  the  following  electrochemical  reactions 
occur  in  the  positive  and  negative  electrodes: 

PbC>2  electrode: 


discharge 

Pb02(s)  +  HSO4  +  3H+  +  2e“  <  »  PbS04(s)  +  2H20 

charge 


Pb  electrode: 


discharge 

Pb(s)  +  HSO4  <  >  PbS04(s)  +  H+  +  2e“ 

charge 

The  governing  equations  with  one-dimensional  assumption 
are  summarized  in  Table  1 .  In  these  equations,  the  effective  prop¬ 
erties,  i.e.  creff,  keff  and  k j^,  are  corrected  to  account  for  electrode 
porosity  as  follows: 

creff  =  <r(l  —  e)ex,  /celT  =  /C£ex,  *{f  =  fcDeex  (1) 

where  (1  —  s)  is  the  volume  fraction  of  conducting  solid  matrix. 

Details  of  the  governing  equations  can  be  found  in  [7]  except 
that  the  equilibrium  potential  Af/pbcb  (at  25  °C)  is  also  added 
to  the  model  from  an  empirical  equation  presented  by  Bode  [8]. 

A£/pbo2  =  1-9228  +  0.147519  log(w)  +  0.063552  log2(m) 

+  0.073772  log3(m)  +  0.033612  log4(m)  (2) 

where  m  is  the  molality  of  the  sulfuric  acid.  Another  empiri¬ 
cal  equation  based  on  literature  data  at  25  °C  is  used  to  relate 
concentration,  c,  to  molality,  m 

m  =  1.00322  x  103c  +  3.55  x  104c2  +  2.17  x  106c3 

+2.06  x  108c4  (3) 


1.2.  Initial/boundary  conditions 

To  solve  the  system  of  equations,  initial  and  boundary  con¬ 
ditions  for  the  primary  variables  are  necessary.  The  initial  con¬ 
dition  for  acid  concentration  is  c  =  c0.  Two  approaches  can  be 
taken  to  find  the  initial  conditions  for  potential  in  solid  and  liquid 
which  are: 

(1)  Solve  the  first  two  steady  equations  given  in  Table  1  with 
constant  c  =  c0. 

(2)  Solve  the  whole  system  with  a  very  small  time  step,  i.e. 
10-4  s. 


Fig.  1.  Schematic  illustration  of  a  lead-acid  cell. 
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The  appropriate  boundary  conditions  at  v  =  0  and  /  for  c  and  0i 
are  ^  ^  =  0;  the  boundary  condition  for  potential  in  solid  is 

0s  =  0  or  V  for  a  given  voltage,  or  —  creff  - ^  =  /  for  a  prescribed 
current  density.  Positive  values  of  I  refer  to  charging,  whereas 
negative  values  of  I  denote  discharging. 


2.  Numerical  procedure 


To  solve  the  system  of  equations  by  the  Keller-Box  method, 
it  should  be  first  converted  into  a  system  of  first-order  partial 
differential  equations  [3].  The  details  of  the  Keller-Box  method 
can  be  found  in  standard  CFD  books  and  have  been  applied 
frequently  in  boundary  layer  theory  [4].  By  defining  ^  =  w, 

^  =  v  and  =  w  and  substituting  into  the  system  of  equa¬ 
tions,  one  can  obtain: 


< 


d(aQifu) 


dx 


=  +Aj 


=  u 


dx 
3(keffn) 
dx 


=  ~Aj  - 


d 

dx 


=  v 


90i 

dx 
dc 

£ —  = 
dt 

dc 


d(Deffw) 
dx 


+  (02 


'  dx 


=  w 


where  c,  0S  and  0i,  related  by  the  well-known  Butler- Volmer 
equation: 


J  =  i  o 


and  overpotential  rj  is  defined  as  77  =  0S  —  0j  —  At/pbCb  for 
positive  electrode  and  77  =  0S  —  0i  for  negative  electrode,  where 
Af/pbo2  is  open  circuit  potential.  The  system  of  Eq.  (4)  together 
with  proper  initial  and  boundary  conditions  mentioned  above 
makes  a  complete  set  of  equations.  To  solve  this  system  of  equa¬ 
tions,  it  should  be  linearized  and  iteratively  solved  (for  example, 
using  Newton  iteration). 

The  transport  equations  for  potentials  in  solid  and  liquid  are 
elliptic  partial  differential  equations  (PDE).  Mathematically,  an 
elliptic  PDE  with  Newman  type  of  boundary  conditions  in  all  of 
the  boundaries  of  domain  has  a  unique  solution  if:  (a)  it  satisfies 
the  compatibility  equation  and  (b)  at  least  one  value  is  known 
inside  the  domain. 

Compatibility  equation  in  battery  dynamics  is  interpreted  as 
the  conservation  of  charge.  It  means  that  the  amount  of  charge 
that  enters  the  cell  at  one  electrode  should  leave  the  cell  at  the 
other  electrode.  To  have  a  unique  solution,  one  should  specify  a 
value  for  potential  in  one  point;  for  example,  0S  =  0  at  the  cen¬ 
ter  of  the  positive  electrode.  Then,  <p\  at  the  center  of  positive 
electrode  can  be  obtained  using  compatibility  equation  [5].  All 
the  potentials  then  are  calculated  related  to  this  reference  poten¬ 
tial.  Without  this  reference  point,  a  unique  solution  cannot  be 
obtained. 


Fig.  2.  Distribution  of  acid  concentration  across  the  cell  during  discharge. 


3.  Results 

The  system  of  governing  equations  has  been  solved  using  the 
Keller-Box  method.  To  verify  the  above-mentioned  procedure, 
the  discharge  problem  of  a  lead-acid  battery  has  been  simulated. 
This  sample  has  been  studied  by  Gu  et  al.  [6]  and  reproduced  by 
Gu  et  al.  [7].  All  the  necessary  parameters  are  the  same  as  the 
ones  used  by  Gu  et  al.  [7]. 

In  Fig.  2,  the  variations  of  acid  concentration  in  time  levels 
0,  60  and  105  s  are  shown.  The  figure  indicates  that  the  results 
of  present  simulation  agree  well  with  the  previous  studies  [6,7]. 
As  it  can  be  seen,  when  the  cell  reaches  cut-off  voltage  (i.e. 
t  =  105  s),  the  acid  is  totally  consumed  in  positive  electrode.  But 
in  the  negative  electrode,  the  acid  is  not  totally  consumed  which 
means  the  negative  electrode  is  overdesigned. 

Fig.  3  shows  the  variations  of  charge  across  the  electrode 
at  different  times.  The  results  show  that  during  the  battery 


Fig.  3.  Distribution  of  charge  across  the  cell  during  discharge. 
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Fig.  4.  Voltage  of  the  cell  during  discharge. 


discharge,  the  amount  of  charge  in  both  electrodes  decreases. 
But  at  the  end  of  discharge,  the  electrodes  still  have  a  lot  of 
charge  which  means  these  electrodes  are  not  fully  utilized. 

Fig.  4  shows  the  simulated  voltage  of  the  battery  cell  versus 
time.  The  results  of  the  present  work  match  very  well  with  the 
results  of  the  other  researchers  [6,7]. 

In  the  previous  works  [6,7],  the  open  circuit  potential 
( A  t/pbo2 ) was  taken  to  be  constant.  To  get  a  more  accurate  simu¬ 
lation,  the  effect  of  acid  concentration  dependency  of  At/pbo2  is 
considered  by  introducing  Bode  relation  [8]  into  the  code.  Fig.  5 
shows  the  comparison  between  the  two  results.  As  it  can  be  seen, 
by  considering  the  effect  of  acid  concentration  dependency  of 
the  open  circuit  potential,  the  voltage  of  the  cell  drops  faster  and 
cut-off  voltage  occurs  at  approximately  90  s  instead  of  105  s. 

Fig.  6  shows  the  acid  concentration  profile  in  both  cases.  At 
the  end  of  discharge,  the  acid  is  not  totally  consumed  in  both 


Fig.  5.  Effect  of  modeling  the  open  circuit  potential  on  cell  voltage  using  Bode 
relation. 


Fig.  6.  Effect  of  modeling  the  open  circuit  potential  on  acid  concentration  using 
Bode  relation. 

electrodes  (comparing  to  the  first  simulation)  and  it  means  that 
both  electrodes  are  overdesigned. 

4.  Conclusions 

The  system  of  transport  equations  of  battery  dynamic  has 
been  solved  numerically  using  the  Keller-Box  method.  The 
results  indicate  that  the  Keller-Box  method  is  a  suitable  method 
for  integration  of  electrochemical  transport  equations  both  in 
integrated  and  multi-region  formulation.  The  boundary  con¬ 
ditions  and  interface  conditions  (in  the  case  of  multi-region 
approach)  can  be  implemented  easily  and  require  no  special  rou¬ 
tine.  The  developed  computer  code  is  capable  of  simulating  lead- 
acid  batteries.  Comparisons  were  made  between  the  simulation 
results  and  the  results  of  other  researchers.  Good  agreement  was 
obtained  and  in  addition,  the  effect  of  acid  concentration  depen¬ 
dency  of  open  circuit  voltage  was  also  investigated. 
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